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ABSTRACT 


We  explore  numerically  the  behavior  of  a  method  of 
describing  the  time  dependent  quantum  mechanics  of  a  curve 
crossing  system.  The  two  nuclear  wave  functions  corresponding  to 
the  two  electronic  states  are  each  described  by  a  Gaussian  wave 
packet.  The  packet  describing  the  Incident  state  mimics  the 
initial  wave  function,  and  the  other  packet  is  created  by  the  time 
dependent  Schroedlnger  equation.  They  are  both  propagated  by 
using  a  variational  method.  The  packets  interact  and  we  do  not 
assume  that  they  have  a  small  width.  Exploratory  calculations  are 
made  for  curve  crossing  dynamics  at  low  kinetic  energy  above  the 
barrier  of  the  lowest  adiabatic  state,  for  tunneling,  for  multiple 
crossings,  and  for  a  curve  crossing  system  which  is  strongly 
coupled  to  a  harmonic  bath  whose  motion  is  described  by  a  mean 
trajectory  classical  Latvgevin  method. 
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I.  INTRODUCTION 

I . 1  Description  of  the  problem  and  its  solution. 

Several  recent  papers^”®  have  shown  that  the  Gaussian  wave 

packet  (GWP)  method,  originated  by  Heller  and  his  coworkers®  and 

6“  1 3 

developed  by  him  and  others  ,  can  be  very  useful  in  treating 

12  2 
atom  diffraction,  diffraction,  rotational  excitation  and 

H-  and  Br_  vibrational  excitation  by  collision  with  a  rigid 

lattice,  as  well  as  atom  diffraction  caused  by  collision  with  a 

4 

surface  undergoing  thermal  motion. 

The  experience  accumulated  so  far  shows  that  the  application 

of  the  GWP  method  to  surface  science  problems  is  very  promising. 

Satisfactory  accuracy  can  be  obtained  with  a  relatively  small 

amount  of  computer  time,  and  the  results  are  easy  to  interpret  in 

terms  of  simple,  intuitive,  classlcal-llke  concepts. 

There  is  however  a  group  of  problems  for  which  the 

application  of  the  GWP  method,  in  the  form  practiced  so  far,  meets 

with  conceptual  difficulties.  A  good  example  is  the  adsorption  of 

a  Li  atom  approaching  a  metal  surface  at  thermal  energies.  It  is 

widely  believed  that  this  process  should  be  described  by  a  curve 

crossing  model,  in  which  Li  sticking  is  a  transition  from  a 

neutral  to  an  ionic  state.  Since  the  motion  of  the  Li  atom  is 

nearly  classical  it  seems  profitable  to  try  to  describe  it  by 

using  Gaussian  wave  packets.  A  correct  GWP  description  of  such  a 

process  requires  the  splitting  of  the  Incident  packet  into  two 

packets,  one  representing  the  ion  bound  to  the  surface  and  the 

other  the  neutral  backscattered  into  the  vacuum.  Packet  splitting 

is  also  required  for  describing  sticking  caused  by  tunneling,  by 

excitation  of  internal  modes  of  the  molecule,  by  excitation  of 

modes  in  the  solid  or  by  transfer  of  energy  from  perpendicular  to 

1-13 

parallel  translational  modes.  The  GWP  methods  employed  so  far 
conserve  the  number  of  packets  and  can  describe  the  above 
processes  only  if  they  are  suitably  modified.  In  this  paper  we 
present  such  a  modification. 

To  describe  the  main  ideas  we  consider  the  case  of  a  Li  atom 
approaching  a  metal  surface.  This  is  characterized  by  a  wave 


function  of  the  form 


i|p(x,R:t)  =  G^(R:t)*j(x,R;t)  +  ( R .  t )  (x ,  R ;  t )  ,  (I.l) 

where  ^j^Cx.R)  Is  the  ionic  state  and  ♦2(x,R)  is  the  neutral  one; 

Gj^(R;t)  and  G2(R;t)  are  the  nuclear  wave  functions  associated  with 

the  electronic  states  and  R  and  x  represent  all  the 

nuclear  and  electronic  coordinates,  respectively.  The  GWP  method 

assumes  then  that  G.  and  G.  are  Gaussian  wave  packets  of  the  kind 
6  ^  ^ 

used  by  Heller: 

G^(R;t)  =  exp{(i/h)[a^{t)(R-R^(t))^+P^(t)(R-R^(t)+y^(t)]} 

(1.2) 

If  a  neutral  Li  atom  approaches  the  surface,  we  have  initially 

Gj^(R;t)aO,  while  the  parameters  in  G2(R;t)  are  determined  by  the 

properties  of  the  incident  atom  (i.e.  direction  of  incidence, 

kinetic  energy,  etc.)  The  approach  of  the  Li  atom  to  the  surface 

can  be  described  by  the  usual  GWP  method,  up  to  the  point  where 

the  transition  to  the  ionic  curve  starts  taking  place.  This 

transition  is  equivalent  to  the  birth  of  a  second  packet,  namely 

G  (R;t).  As  the  time  evolves  the  two  packets  interact  with  each 

2 

other  and  build  up  the  correct  ionization  probability  |Gj^{R;t)| 

To  handle  this  situation  the  GWP  procedure  must  be  extended 
in  several  ways.  First,  we  must  find  a  way  of  creating  G^^,  at  the 
proper  moment;  we  use  for  this  a  short  time  Green's  function  which 
generate  G^^  from  G2,  through  the  coupling  between  the  ionic  and 
the  neutral  states.  Once  the  new  packet  is  created,  it  must 
interact  with  the  initial  one,  so  that  their  joint  evolution  leads 
to  the  correct  ionization  probability.  This  is  achieved  by  using 
the  minimum  error  method.^®  Finally,  we  do  not  assume  that  the 
packets  are  narrow  throughout  the  interaction  region  since  this 
approximation  exagerates  the  classical  character  of  the  Li^  motion 

Q  1  i 

and  leads  to  errors .  ' 


The  resulting  MEM  split  Gaussian  wave  packet  (MEM-SGWP) 
method  can  be  used  to  explore  several  physical  and  computational 
issues  pertinent  to  the  above  model.  First  we  investigate  whether 
the  split  packets  behave  reasonably  and  whether  SGWP  is 
numerically  stable;  then  we  examine  several  simplified  versions  of 
the  theory  to  see  whether  they  work  well.  Numerical  studies  meant 
to  answer  such  questions  are  presented  in  Section  III.  2. 

A  more  subtle  and  uncertain  matter  is  whether  SGWP  can 
generate  tunneling  behavior,  that  is,  whether  an  Incident  packet 
whose  energy  is  below  the  height  of  the  barrier  on  the  lowest 
adiabatic  potential,  can  be  split  to  create  a  packet  behind  the 
barrier.  The  customary  GWP  method  will  not  permit  tunneling  since 
the  center  of  the  Incident  packet  moves  on  a  classical  trajectory. 
Calculations  demonstrating  tunneling  behavior  are  presented  in 
Section  III. 3. 

Further  numerical  calculations  (Section  III. 4)  attempt  to 
establish  when  a  two  state  description  is  necessary  and  what  are 
the  consequences  of  this  necessity.  Even  in  those  cases  (e.g. 
charge  transfer)  where  a  two  dlabatic  state  description  is 
convenient  it  does  not  always  follow  that  a  two  state  description 
is  also  needed  in  the  equivalent  adiabatic  representation.  We 
expect  to  need  two  adiabatic  states  when  the  energy  of  the 
incident  packet  is  comparable  with  gap  between  the  states.  In 
such  cases  the  "transmission"  of  the  wave  function,  through  the 
region  where  the  two  adiabatic  states  are  closest  to  each  other, 
is  impaired  as  if  the  upper  state  helps  "reflect"  the  incident 
packet.  Calculations  showing  that  SGWP  generates  such  behavior 
are  presented  in  Section  III. 4. 

In  a  one  dimensional  system  the  creation  of  a  packet  on  the 
ionic  state  does  not  lead  to  binding:  the  ionic  packet  oscillates 
in  the  well  and  splits  into  an  outgoing  neutral  packet  and  an 
ionic  packet  whenever  it  goes  through  the  crossing  region. 

Repeated  occurrence  of  this  process  makes  the  amplitude  of  the 
nuclear  ionic  state  go  to  zero.  The  above  sequence  is  the  GWP 
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description  of  a  predissociation  process.  The  calculations 
presented  in  Section  IV  show  that  a  multiple  crossing  SGWP  model 
leads  to  such  behavior. 

Finally  in  Section  V  we  explore  the  curve  crossing  behavior 

of  a  neutral  atom  which  is  ionized  upon  approaching  a  moving 

lattice.  The  manner  in  which  the  thermal  lattice  motion  modifies 

the  behavior  of  the  packets  is  described  by  coupling  the  MEM-SGWP 

equations  to  a  mean  trajectory  classical  Langevin  equation 

describing  lattice  motion.  In  surface  science  such  a  model  is  . 

relevant  with  regard  to  sticking  of  alkali,  which  is  believed  to 

require  a  two  state  description;  previous  work  on  adsorption- 

14  15 

desorption  dynamics  '  considered  one  energy  surface  only.  The 
calculation  is  also  relevant  to  the  problem  of  tunneling  in 
systems  subject  to  thermal  noise,  which  has  received  a  lot  of 
attention  lately.^® 

All  the  calculations  presented  here  are  exploratory  and 
intend  to  test  the  qualitative  behavior  of  SGWP.  We  plan  to  study 
the  accuracy  of  the  method  once  we  develop  an  exact  procedure  for 
solving  the  problems  described  above. 

I . 3  Other  methods 

Most  of  the  existing  work  applying  curve  crossing  models  to 

17-32 

surface  science  problems  has  often  been  concerned  with  high 

kinetic  energy  phenomena,  where  sticking  is  not  an  issue  and 
classical  approximations  of  varying  quality  are  fairly  adequate. 
Here  we  are  especially  concerned  with  the  limit  of  very  low 
kinetic  energy  where  single  classical  trajectory  methods  are  not 
likely  to  be  useful  and  quantum  effects  (both  for  the  electronic 
excitation  and  perhaps  the  motion  of  the  light  atom)  are 
important . 

There  are  a  number  of  approaches,  other  than  GWP ,  that  might 

be  usefully  applied  to  the  problem  of  interest  here.  The  mean 

28-29  32-39 

trajectory  approximation  (MTA)  '  uses  one  classical 

trajectory  to  describe  the  motion  of  the  nucleus.  The  two  state 
aspects  of  the  problem  are  partly  incorporated  by  using  in  the 
classical  equation  an  average  potential  energy  which  includes 


contributions  from  both  the  ionic  and  the  neutral  surfaces,  in 

proportion  to  their  instantaneous  occupation  (see  Section  III  for 

39 

details).  Recent  calculations  show  that  at  low  kinetic  energies 

the  method  leads  to  unphysical  behavior. 

WKB  theory  has  been  extended  to  curve  crossing 
40-43 

problems  and  its  main  limitation,  as  far  as  surface  science 

is  concerned,  is  the  difficulty  of  extending  it  to  three 
dimensions . 

44 

The  multiple  trajectory  method  of  Tully  and  Preston  is 

frequently  critized  because  it  neglects  interference  effects; 

since  these  are  likely  to  be  washed  out  by  thermal  averaging  (when 

coupling  to  phonons  is  included)  this  is  not  a  significant  flaw 

here.  A  potentially  important  limitation  is  that  it  will  miss 

quantum  effects  for  light  particles. 

Another  fruitful  line  of  research  was  originated  by  attempts 

to  Implement  and/or  simplify  the  semi-classical  method  developed 

45 

in  several  elegant  papers  of  Pechukas.  Straight  applications  of 

Pechukas'  method  led  to  inefficient  numerical  codes  and  numerical 
46 

instabilities.  The  simplifications  introduced  by  George  and 
47 

Miller  offer  some  remarkable  insights  in  the  physics  of  the 

problem,  but  have  not  yet  reached  the  computational  simplicity 

characterizing  the  GWP  approach. 

Finally  Herman  and  Freed's  implementation^®^  of  early  work 
49b 

of  Laing  and  Freed  leads  to  an  interesting  method  which  however 
seems  rather  difficult  to  apply,  especially  for  many  dimensions. 

In  this  context  SGWP  has  several  appealing  features;  it  .is 
easy  to  use  in  three  dimensions;  it  is  computationally  efficient: 
and  it  provides  an  Intuitive  classical-like  description  of  time 
dependent  quantum  mechanical  process. 
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II.  THE  MODEL 

II . 1  The  Hamiltonian 

49 

We  consider  a  system  described  by  two  diabetic  electronic 
states  <l>^(x,R)  ,  i  =  1,2  and  the  Hamiltonian  matrix 


H^j(R)  =  J-dx*  .  (x,R)H{x.R)<I>j(x,R) 


{II. 1) 


with 


lii(R)  =  W^{exp[-2a(R-R^) ]  +  (-1)^  exp[ -a( R-R^ ) ] } 
+  (-l)^~^Ae/2 


II  .  2  ' 


and 

H^j  =  (pAe/2)  exp{-o^(R-Re)^}  (II. 3) 

The  "ionic  curve"  is  a  Morse  potential  whose  asymptotic  energy 

(i.e.  the  energy  for  R-*<»)  is  not  zero,  as  customary,  but  Ae/2. 

The  asymptotic  energy  gap  between  the  two  states  (i.e.  the 

"ionization  potential")  is  Ae.  The  binding  energy  of  the  ion  to 

the  surface  is  and  p  is  a  constant  characterizing  the  strength 

of  the  coupling  between  the  states.  The  off  diagonal  element 

peaks  at  R=R  which  is  the  position  where  H, ,  and  H, _  cross. 

The  values  of  the  parameters  used  in  the  calculations 

presented  here  are  listed  in  Table  I.  The  corresponding 

Hamiltonians  are  called  Hamiltonian  I  and  II.  In  Figs.  1,  and  4 

we  plot  the  diabatic  potential  surfaces,  the  adiabatic  ones  and 

the  coupling  H^2  these  Hamiltonians. 

Since  i  =  l,2  are  diabatic  states  the  coupling  between 

them  takes  place  through  H  •  therefore,  for  simplicity,  we  ignore 

^  ^49 

other  coupling  terms  and  set 


J*dX  (  00  ,  R) 


<I»j  (x,R) 


6  .  .6 
ij  no 


(II. 4) 


Using  the  wave  function  (I.l)  and  the  conditions  (II.  4)  in 
the  time  dependent  Schrodinger  equation  leads  to 
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ih3G^(R,t)/3t  =  (K  +  H..(R))G^  +  H.^G^  ,  (II. 5) 

where  K  is  the  kinetic  energy  operator.  We  solve  this  equation 

under  the  assumption  (which  is  the  essence  of  Heller's  method) 

that  the  functions  G.  have  the  form  (1.2)  and  maintain  it  at  all 

1 

times;  the  effect  of  the  collision  is  to  change  the  parameters 
appering  in  (1.2).  In  the  present  work  we  derive  differential 
equations  and  initial  conditions  for  ot^,  R^,  and  and  solve 
them  numerically  for  the  problems  specified  in  Section  I. 

II . 2  The  initial  conditions 

II. 2. a.  The  creation  of  a  new  packet 

Consider,  for  illustration,  the  case  when  the  colliding 
particles  are  initially  in  the  electronic  state  2.  The  parameters 
of  the  packet  G2(R:t)  are  determined  by  the  Initial  conditions, 
which  also  set  Gj^(R;t)  =  0. 

Since  all  Heller-like  methods  conserve  the  number  of 
Gaussians,  we  must  make  a  modification  that  will  split  the 
incident  packet  into  a  "neutral"  and  an  "ionic"  one.  To  achieve 
such  splitting  we  use  the  formal  solution  of  Eq .  (II. 5): 

|G^(t)>=|G^{t^)>-(i/h)  X  dt  exp{-(i/h)[K+H^^J(t-t' ) 

^o 

(II. 6) 

Here  t^  is  the  time  at  which  the  packet  begins  to  overlap 

spatially  with  This  is  the  last  time  step  in  the  integration 

of  the  differential  equation  (II.  5)  for  which  we  can  maintain 
|Gj(t^)>=0.  In  the  next  time  step  we  must  create  a  new  packet 
|G^(t^-t-T)>  whose  functional  form  is  determined  from  (II.  6).  By 
using  a  very  small  value  for  t  we  can  replace  the  integral  with  x 
exp { - ( i /h ) [ K+H  ] T }H  I G  ( t  )>  and  the  exponential  operator  with 
its  short  time  limit.  This  gives 

Gi(R;to+T)  =  -(  iT/h)  (m/2TTihT  )  ^-^^ 


fc  "  1.  *  ^ 

* 


XdR'exp{ (i/h) [ (m/2T) (R-R- ) ] } 


H2j{R' )Gj(R;t^) 


(II. 7) 


Generally  the  function  Gj(R;t^+T)  is  not  a  Gaussian. 

However,  since  t  Is  arbitrarily  small  we  can  use  -  with  arbitrary 
accuracy-  the  stationary  phase  approximation  which  replaces 
Hjj(R')  with  it's  second  order  expansion  in  powers  of  (R-R' ) . 

Since  in  the  present  model  H2j^(R')  and  Gj^(R';t^)  are  both 
Gaussians,  the  use  of  the  above  expansion  allows  us  to  perform  the 
integral  in  (11.7)  analytically;  this  gives  for  Gj^(R;t^+T)  a 
Gaussian. 

We  can  summarize  this  whole  procedure  by  stating  that  we 

generate  G^^  from  G^  by  using  first  order  perturbation  theory  with 

respect  to  the  small  parameter  tH..,  and  a  stationary  phase 

2 

approximation  with  respect  to  the  large  parameter  ml  /hi,  where  1 
is  a  length  over  which  the  potential  changes.  Both 

approximations  are  arbitrarily  accurate  since  t  is  arbitrarily 
small . 

II. 2. b  The  initial  parameters  of  the  incident  packet 

Since  in  all  GWP  methods  the  parameters  of  the  Gaussian  wave 
packets  evolve  according  to  first  order  differential  equations, 
the  theory  is  not  defined  unless  it  provides  a  prescription  for 


the  choice  of  the  initial  conditions. 


In  the  GWP  method  this 


simple  requirement  can  become  a  delicate  and  ambiguous  matter. 

In  general,  the  time  dependent  theory  of  collisions  uses  an 
initial  wave  function  which  mimics  the  pre-collision  state 
prepared  by  the  experimental  set  up.  This  causes  difficulties  for 
the  GWP  theories  since  it  is  unlikely  that  anyone  would  want  to  do 
-  at  least  in  the  foreseable  future  -  experiments  in  which  the 
initial  state  is  a  narrow  wave  packet.  To  understand  why  this  is 
the  case  we  consider  how  an  idealized  experiment  of  this  kind 
might  be  done.  We  can  prepare  a  wave  packet  by  measuring  the 
position  of  the  "projectiles"  (i.e.  the  incident  molecules)  before 
they  reach  the  surface.  To  do  this  we  must  intersect  the 
projectiles  beam  with  a  probe  beam  designed  to  overlap  with  it  in 


a  very  small  volume  1  .  In  such  an  arrangement  the  deflection  of 
a  probe  particle  at  the  time  t^  signals  the  fact  that  a  projectile 
was  present  in  the  volume  1^.  The  state  of  that  projectile  can  be 
adequately  described  by  a  packet  of  width  1,  whose  center  was 
located  in  the  probe  volume  at  t^ .  The  packet  thus  formed 
continues  to  travel  towards  the  surface,  collides  with  it  and  is 
scattered  into  the  detector;  the  particles  in  a  plane  state  - 
which  went  undisturbed  through  the  volume  -  have  the  same  fate. 

The  detector  registers  the  arrival  of  both  the  packet  and  the 
plane  wave  and  cannot  distinguish  between  them.  To  measure  only 
the  arrival  of  the  deflected  wave  packets  we  must  lower  the  flux 
of  incident  projectiles  to  the  point  that  whenever  a  projectile 
scattered  by  the  surface  is  detected  in  coincidence  with  the 
deflection  of  a  probe  particle  we  can  be  fairly  certain  that  we 
are  dealing  with  packet  scattering. 

Since  there  is  a  very  low  probability  that  such  difficult 

experiments  will  be  performed  we  must  decide  what  is  the 

experimental  significance  of  these  GWP  calculations.  One  point  of 

view^^  is  that  in  all  experiments  -  other  than  the  one  described 

above  -  the  packets  have  no  reality.  They  are  merely  "pieces"  of 

the  wave  function,  Introduced  for  computational  convenience,  and 

only  the  coherent  sum  of  such  packets  -  which  represents  the  total 

wave  function  -  has  any  meaning.  Therefore  the  choice  of  the 

initial  parameters  in  each  Gaussian  G  should  be  made  so  that  Z  G 

l8  “ 

best  fits  the  initial  wave  function. 

Another  point  of  view  is  that  a  Gaussian  wave  packet  state 
is  a  limiting  case  in  quantum  mechanics,  which  we  might  call 
corpuscular  (as  opposed  to  wave  like),  which  provides  a  reasonable 
description  of  a  semi-classical  time  dependent  processes.  This  is 
the  point  of  view  which,  for  the  sake  of  simplicity,  is  adopted 
here . 

Unfortunately  this  interpretation  does  not  provide  a 
prescription  for  choosing  the  initial  values  of  all  packet's 
parameters.  We  can  choose  R  (t=0)  so  that  the  packet  is  just 
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outside  the  interaction  zone;  we  can  select  P2(t=0)  /2m  to  give 

the  initial  energy  (there  is  a  slight  error  involved  in  doing 

this) ;  we  can  choose  a  normalized  packet  and  thus  have 

exp{-2Imy2/^}  (J'fT/2Ima2)  =  1.  Since  the  initial  phase  of  the 

packet  is  irrelevant  we  can  choose  Rey2(^“0)“0'  We  are  thus  left 

with  two  unknown  parameters,  Rea2(t=0)  and  Ima2(t=0), 

To  get  the  calculation  started  one  customarily  follows 
0 

Heller  who  recommends  a  choice  of  parameters  that  would  give  a 
narrow  minimum  uncertainty  wave  packet  at  the  most  important  point 
on  the  trajectory.  In  potential  scattering  this  point  is  the 
turning  point.®  In  the  present  case  there  is  some  ambiguity  since 
we  can  choose  either  the  crossing  point  or  the  turning  point, 
which  are, both  Important.  We  have  selected  the  former  on  the 
subjective  belief  that  the  transfer  from  one  curve  to  another  is 
the  most  important  event  in  the  present  system.  Therefore  we  use 
the  equations  of  motions  for  Rea(t)  and  Ima(t)  to  select  the 
initial  values  Rea(t=0)  and  Ima2(t=0)  that  will  give  Reot2(t^)=0 
and  a  small  width  l-(t  )  =  (h/2Imoc-(t  at  the  time  t  at 

C  2  C  C 

which  the  center  of  the  packet  G2  reaches  the  crossing  position 

R^ .  Specifically  we  chose  Rea2(t^)=0  and  a  large  value  for 

Ima  ( t  )  and  perform  a  backward  propagation  from  R  to  the  initial 
&  ^  c 

position  R2(0)  ,  by  using  the  equations  of  motion  for  Rea2  In'«2 

on  the  surface  H22(R)<  without  any  coupling  to  the  state  1. 

We  note  that  we  see  no  compelling  a  priori  reason  to  follow 

Heller's  recommendations.  There  is  some  a  posteriori 

justification  since  in  the  past^  ®  such  choices  led  to  accurate 

results  which  were  stable  with  respect  to  variations  of  the 

1-4  6 

initial  values  chosen  for  0!(t=0).  ’  We  find  that  the  same 

stability  exists  in  the  present  system  but  we  do  not  necessarily 
expect  this  to  be  the  case  for  new  problems  of  a  different 
character.  At  this  time  we  regard  this  ambiguity  as  a  temporary 
nuisance  which  can  be  avoided  by  using  GWP-s  to  fit  the  initial 
wave  function  that  is  prepared  experimentally. 

II . 3  The  propagation  of  the  wave  packets 

Once  the  parameters  of  the  initial  G2  packet  are  chosen  and 


r  j 
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the  new  G.  packet  is  created  we  must  select  a  scheme  for  their 
propagation.  The  simplest  Heller  scheme  assumes  that  the  packets 

remain  narrow  throughout  the  collision,  and  move  independently. 

8  9 

Skodje  and  Truhlar,  and  Heather  and  Metiu  have  shown  that  the 
first  assumption  does  not  work  well  for  packets  which,  like  , 
are  trapped  in  an  anharmonic  well.  The  second  assumption  -  which 
is  reasonable  for  the  problems  that  were  of  interest  to  Heller®  - 
is  untenable  in  the  present  case,  since  the  amplitude  of  the 
packet  Gj^  must  grow  at  the  expense  of  that  of  ;  this  cannot  be 
achieved  if  the  packets  are  decoupled. 

For  these  reasons  we  are  using  for  the  propagation  a 
variational  procedure which  we  call  the  minimum  error  method 
(MEM)  which  does  not  make  the  assumptions  mentioned  above.  The 
only  approximation  is  that  in  the  course  of  time  the  wave 
functions  G^  and  G2  remain  Gaussians. 

In  working  with  other  systems  Sawada,  Heather,  Jackson  and 
Metiu^^  have  shown  that  sometimes  the  use  of  one  Gaussian  per 
electronic  state  does  not  offer  enough  flexibility  in  the 
variational  wave  function;  improved  accuracy  can  be  obtained  by 
using  nuclear  wave  functions  which  are  sums  of  Gaussians  coupled 
to  each  other.  We  have  developed  such  a  method  for  the  curve 
crossing  problem  and  preliminary  numerical  tests  show  it  to  be 
rather  expensive.  Therefore  at  this  time  we  prefer  to  use  the 
present  one-Gaussian-per-electronic-state  model . 

The  variational  equations  obtained  by  using  the  minimum 
error  method  are  given  in  the  Appendix.  The  equations  are 
transformed  by  using  Heller's  P-Z  method  and  the  resulting 
differential  equations  are  solved  by  using  a  fourth  order  Runge- 
Kutta  method  with  variable  step  size. 

I I . 4 .  Numerical  details  and  the  units. 

All  the  calculations  presented  here  were  carried  out  with 
the  Hamiltonians  defined  by  Eqs .  (II. 1-3)  and  the  parameters  given 
in  Table  I.  We  use  a  "natural"  system  of  units  with  ot  ^  for 
length,  («v^)  for  time  and  hocv^”^  for  mass;  the  parameter  ot”^ 
controls  the  range  of  the  potentials  (see  Eqs.  II.  (2-3))  and  v 
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is  the  initial  velocity.  The  latter  is  the  expectation  value  of 

the  velocity  operator  -(ih/M)3/3R  for  the  initial  wave  packet  and 

it  is  equal  to  P2{t)/M,  where  P2(t)  is  the  momentum  of  the  packet 

(see  Eq.  1.2).  An  important  quantity  is  Massey's  parameter  X  = 

Ae/hotv^,  which  is  the  asymptotic  energy  gap  Ae  (i.e.  the 

ionization  potential)  in  the  energy  units 

The  matrix  elements  corresponding  to  the  parameters 

given  in  Table  I  are  shown  in  Pigs.  1.  4  and  8,  together  with  the 

corresponding  adiabatic  states.  On  the  same  figures  we  have 

indicated  by  arrows  the  kinetic  energies  of  the  incident  packets 

used  in  the  present  calculations.  These  are  labelled  by  the 

values  of  the  Massey  parameter  X,  whose  inverse  is  proportional  to 

the  incident  velocity.  Note  that  for  a  particle  in  a  GWP  state  G^ 

the  kinetic  energy  <G- | ( -h^/2m) 3^/3R^ | G*>  differs  from  the 

classical  energy  P-(t)  /2m  of  the  center  of  the  packet.  The 

^  2 

difference,  denoted  AP-/2m,  is  shown  in  Table  II,  which  also  shows 

2  ^ 

the  values  of  X,  of  p2/2m  and  the  difference  between  the  kinetic 

energy  of  the  Incident  packet  and  the  height  of  the  barrier  on  the 

2 

lowest  adiabatic  state.  Note  that  AP-/2m  is  very  small  compared 

^  2 

to  the  "classical"  kinetic  energy  of  the  packet  P2(t=0)  /2m, 
because  we  use  a  spatially  broad  packet  which  has  a  narrow  spread 
in  the  momentum  space. 


Ill .  PACKET  SPLITTING  FOR  ONE  CROSSING 

III.l  The  behavior  of  SGWP  method  and  Its  comparison  with  MTA  and 
LHA. 

In  this  section  we  report  calculations  with  the  Hamiltonian 

I,  defined  by  Eqs .  (11.1-3)  and  the  parameters  listed  in  Table  I. 

The  matrix  elements  as  well  as  the  corresponding  adiabatic 

states  are  shown  in  Fig.  1.  We  examine  the  behavior  of  the 

packets  generated  by  the  split  GWP  method  and  compare  them  to  that 

given  by  simplified  versions  of  the  theory. 

Since  we  compare  the  present  calculations  (MEM)  with  the 

mean  trajectory  approximation  (MTA)  and  the  local  harmonic 

approximation  (LHA)  we  describe  them  briefly  here.  MTA  assumes 

2  8 

the  wave  function 


4((x,R;t)  =  G(R;t)  [c^(t)<t^(X.R)  +  C^  ( t )  *2  ( X ,  R )  ]  ,  (III.l) 

2  8 

where  G  is  a  GWP  of  the  form  (1.2).  One  can  show  that  the 

center  of  the  Gaussian  G  moves  according  to  Newton’s  equation  with 

« 

the  force  -S  Z  c^ ( t ) c ^ ( t ) 3H^ ^ ( R) /3R.  While  this  force  depends  on 
the  Instantaneous  populations  of  the  two  states  and  the  phase 
difference  between  the  corresponding  amplitudes,  it  is  unable  to 
generate  a  nuclear  motion  which  has  a  "two  trajectories" 
character.  This  shortcoming  is  important  only  at  low  incident 
kinetic  energy. 

The  local  harmonic  approximation  (LHA)  uses  one  packet  for 
each  electronic  state  and  assumes  that  G^^  and  G2  are  narrow 
throughout  the  collision  so  that  the  integrals  in  which  they  are 
involved  can  be  computed  by  the  method  of  steepest  descent  or  the 
stationary  phase  approximation.  This  simplifies  the  differential 
equations  which  give  the  evolution  of  the  parameters  in  G^^  and  G2 
and  speeds  up  their  integration.  Physically  the  use  of  narrow 
Gaussian  amounts  to  taking  a  semiclassical  limit. 

In  all  the  calculations  presented  here  we  assume  that  the 
initial  state  is  neutral .  The  packet  G2  is  constructed  as 
discussed  in  Section  II. 2b  and  the  packet  G,  is  created  as 


described  in  Section  II. 2a. 

In  Fig  2  we  show  how  the  centers  of  the  two  packets  move, 
for  various  incident  kinetic  energies  of  .  The  time  evolution 
of  the  corresponding  state  occupation  probabilities  is  shown  in 
Figs .  3 . 

Let  us  examine  in  detail  Fig.  2a  which  Illustrates  the 
motion  of  the  two  packets.  At  the  initial  time  the  center  of  the 
packet  G2  is  located  at  =  10,  and  that  of  the  newly  created  G^^ 
is  at  Rj  =  8.2.  The  new  packet  is  closer  to  the  crossing  point 
located  at  R  =5  (see  Fig.  1))  than  the  incident  one,  but  they  are 
both  far  from  it.  After  about  5  time  units  the  packet  G^  catches 
up  with  G^ •  In  all  the  calculations  carried  out  so  far  the  new 

packet  G.  is  created  closer  to  the  crossing  point  and  G.  catches 

^  2  2^ 
up  with  it  because  its  energy  of  (i.e.  <G2l-h  ^  /2m|G2>)  at  the 

moment  of  G^^'s  birth  is  larger  than  that  of  G^^ .  At  the  time  6.5 

the  packet  G2  reaches  its  turning  point,  i.e.  its  center  turns 

around  and  moves  away  from  the  surface.  Note  that  the  turning 

point  of  the  center  (at  R=5.4)  is  not  the  classical  turning  point. 

This  behavior  is  different  from  that  given  by  the  simplest  version 

of  Heller's  theory  -  which  makes  the  local  harmonic  approximation 

(LHA)  -  in  which  the  center  of  the  packet  moves  classically  and 

turns  at  the  classical  turning  point.  As  shown  by  Heather  and 

Metiu'*’  if  LHA  is  not  made  the  center  of  the  packet  does  not  move 

classically  and  behaves  like  a  "fuzzy  ball"  which  turns  upon 

collision  with  a  repulsive  potential  before  its  center  reaches  the 

wall.  Furthermore,  the  coupling  between  Gaussians  introduces  a 

new  force  in  the  equation  of  motion  for  R2(t),  which  has  no 

classical  analog.  This  can  also  affect  the  trajectory  and 

1 2 

the  location  of  the  turning  point. 

While  the  packet  G2  turns  around  and  leaves  the  scene  of  the 
action,  the  center  of  G^^  is  accelerated  towards  the  wall  of  the 
potential  turns  around  at  R^^  2:  1.2,  is  trapped  in  the  well 

and  oscillates  in  it.  The  motion  of  this  packet  does  not  have  to 
be  classical  but  it  happens  to  resemble  classical  motion.  In 
the  cases  in  which  the  energy  of  G  exceeds  the  dissociation 


energy  of  the  packet  escapes  from  the  surface. 

The  Figs.  2b-d  correspond  to  different  kinetic  energies  of 
the  incident  packet  and  show  the  same  general  behavior:  the  split 
packets  move  in  a  nearly  classical  fashion  on  the  two  surfaces. 

By  comparison  the  mean  trajectory  approximation  fairs 
poorly.  The  packet  G  representing  the  incoming  particle  manages 
to  penetrate  above  the  ionic  well  and  for  a  brief  period  it  moves 
similarly  to  G^.  However  when  it  reaches  the  crossing  point  for 
the  second  time  it  changes  its  behavior  to  resemble  G^^.  This 
happens  because  at  t=5  the  effective  potential  =2!  Z  c^  ( t )  * 

Cj(t)H^j(R)  switches  from  resembling  to  resembling 

(Ic^l^  becomes  small  and  grows,  as  seen  in  Fig.  3a.)  At  the 

time  t=ll  the  behavior  is  reversed  because  IC2I  starts  growing 
rapidly  again  and  ~  H22* 

Since  packet  trajectories  are  not  observables  we  should 

judge  the  usefulness  of  MTA  by  comparing  its  predictions  for  the 

occupation  probabilities  to  those  of  MEM-SGWP.  This  is  done  in 

Figs.  3a-d.  The  results  obtained  with  MTA  are  consistently  poor. 

We  note  that  a  detailed  numerical  analysis  of  MTA  shows  that  at 

3  3 

low  kinetic  energies  MTA  has  rather  unphysical  behavior. 

The  LHA  approximation  for  the  split  Gaussians  does  give 
trajectories  which  are,  qualitatively  speaking,  well  behaved. 
However  the  asymptotic  values  of  the  transition  probability  are 
different  from  those  of  MEM-SGWP  as  shown  in  Table  III. 

We  have  found  that  all  the  results  reported  here  are  stable 
with  respect  to  the  choice  of  the  initial  width  of  G^  and  of  the 
point  where  G^  is  first  created.  If  G^  is  created  too  early  (i.e, 
at  a  distance  which  is  too  far  from  the  point  where  starts 

being  different  from  zero)  subsequent  propagation  makes  it 
disappear  (i.e.  its  amplitude  goes  to  zero).  If  it  is  started  too 
late  (i.e.  in  the  crossing  region)  the  results  are  strongly 
dependent  on  the  starting  point.  If  G^  is  created  in  a  strip  of 
about  1 . 5A ,  located  far  from  the  crossing  point  but  in  a  region 
where  H^^  is  not  exactly  zero,  its  properties  are  independent  on 
the  point  of  creation.  In  all  cases  we  find  that  immediately 


after  the  appearance  the  parameters  in  the  Gaussians  undergo 

fast  transient  changes  and  then  settle  to  a  smooth  and  physically 
reasonable  evolution. 

Ill, 3  Tunneling  within  SGWP 

In  this  section  we  study  systematically  how  the  packets 
behave  as  the  energy  of  the  incident  packet  changes  from  being 
above  to  being  below  the  barrier  on  the  lowest  adiabatic  state. 

The  calculations  are  carried  out  with  the  Hamiltonian  II  defined 
by  Eqs.(II.l-3)  and  the  parameters  given  in  Table  I,  The  value  of 
^  is  0.3  and  the  mass  is  that  of  Li.  The  matrix  elements  H^.{R) 
are  plotted  in  Fig,  4  together  with  the  potential  energies  of  the 
corresponding  adiabatic  states.  The  incident  kinetic  energies  are 
labelled  by  the  values  of  X  and  are  indicated  on  the  graph.  The 
difference  between  the  incident  kinetic  energy  and  the  top  of  the 
barrier  on  the  lowest  adiabatic  state  is  given  in  Table  II. 

The  trajectories  followed  by  the  MEM-SGWP  packets  is  shown 
in  Figs.  5a-f.  The  general  behavior  parallels  that  already  seen 
in  the  sequence  shown  in  Fig.  2.  In  the  first  three  calculations 
(i.e.  X  =  100,  110,  120)  the  kinetic  energy  is  above  the  barrier 
and  in  the  last  three  (Figs.  (6d-f))  it  is  below.  The 
trajectories  are  only  slightly  modified  as  the  incident  energy  is 
slightly  lowered,  with  no  apparent  trauma  when  the  energy  gets 
below  the  barrier.  The  lowest  energy  curve  is  somewhat  peculiar 
since  both  packets  turn  around.  This  is  not  acceptable  behavior 
for  the  ionic  packet  G^^,  but  it  is  not  cause  for  concern:  the 
amplitude  of  the  misbehaving  packet  is  completely  negligible! 

It  is  more  interesting,  from  the  point  of  view  of  tunneling, 
to  monitor  the  total  energies  of  each  packet  as  a  function  of 
time.  At  this  point  it  is  necessary  to  make  a  few  remarks 
concerning  the  energy  in  the  SGWP  theory. 

The  total  energy  of  the  system  is 

H(t)=N(t)  ;  dxdR  ^^l•(x,R;t)H  i)i(X,R;  t)  (III.l) 


with  the  normalization 
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N(t)  H  ;dxJ‘dRil)*(x,R;t)i|)(x,R:t)  =  Z  /dRG  .  (R;  t )  G  (R;  t )  .  (III. 2) 

1  ^  ^ 

By  using  (I.l)  in  Eq.  (III.l)  we  can  write 

H(t)=N(t)"^Z  J‘dRG*{K+H^^(R)  } G^+N ( t )  “^2ReJ‘dRG*  ( R ;  t )  2  < ^ ^ 

(III. 3) 

The  terms  in  the  sum  represent  the  results  that  would  be  obtained 
in  a  state  selected  energy  measurement  (e.g.  Ej^(t)  is  obtained  if 
the  energy  of  the  ionic  component  is  measured) .  The  last  term  in 
(III. 3)  is  an  energy  contribution  due  to  quantum  interference 
between  the  nuclear  wave  functions. 

In  Figure  6  we  plot  the  time  evolution  of  the  energies 


E^(t)  =  /dRG^[K+H^^(R) ]G.//dRG^G^ 


(III. 4) 


which  differ  from  the  terms  appearing  in  the  sum  in  Eq.  (III. 3) 
only  through  normalization.  We  believe  that  plotting  E^(t)  as  a 
function  of  time  gives  a  feeling  for  the  rate  of  energy  exchange 
between  the  packets;  furthermore  the  position  of  Ej^(t)  and  E2(t) 
with  respect  to  the  barrier  on  the  lower  adiabatic  states  can  be 
used  as  an  indication  whether  tunneling  takes  place. 

In  Fig. 6a  we  show  the  evolution  of  E^(t)  for  X=100.  This 
corresponds  to  a  kinetic  energy  above  the  barrier  on  the  lowest 
adiabatic  state  and  below  the  energy  of  the  upper  adiabatic  state 
(see  Fig.  5).  The  energy  of  the  newly  formed  packet  at  the 
time  of  formation  (i.e.  t=7  is  much  higher  than  that  of  the 
incident  packet.  This  does  not  cause  any  problem  with  energy 
conservation.  To  show  this  we  write  the  total  energy  as 


E(t)  =  Z  IP.(t)E.(t)  +  N(t)'^2ReXdRG*H^2^2 
where 

*  2 

(t)  =  XdRG  G  /  Z  XdRG.G 
^  ^  1=1  ^  ^ 


(111.5) 


III .6) 


is  the  probability  that  the  particle  is  in  the  state  i.  The  total 


energy  E(t)  is  conserved  at  all  times  even  when  Ej(t)  is  very 
large,  because  when  Ej{t)  grows,  ^^(t)  becomes  small;  furthermore 
the  interference  term  can  also  compensate  some  of  the  changes  in 
Ej^{t)  to  keep  E(t)  constant.  The  same  kind  of  reasoning  explains 
how  it  is  possible  that  total  energy  is  conserved  in  the 
asymptotic  channels  (i.e.  bound  and  unbound  G^)  even  though  the 
packet  energies  E^  are  changed  by  the  collision. 

We  note  several  interesting  results  concerning  the  energy  of 
the  packets.  In  all  our  calculations  the  energy  of  the  newly 
formed  packet  is  above  the  asymptotic  value  H^^(R-«o).  Thus  the 
center  of  the  packet  is  not  placed,  when  it  is  created,  in  a 
classically  forbidden  region  of  the  diabatic  state  In 

fact,  in  all  cases  in  which  the  packet  G^^  manages  to  penetrate  in 
the  region  above  the  ionic  well  it  does  so  without  going  through  a 
classically  forbidden  region.  We  emphasize  that  in  our 
propagation  scheme  we  do  not  use  a  quadratic  expansion  of  the 
potential  and  do  not  assume  independent  Gaussians.  Thus  neither 
the  motion  of  the  center  of  the  packet  (i.e.  the  evolution  of 
R^(t)  and  P^(t))  nor  the  energy  of  the  packet  are  classical,  so 
there  is  no  a  priori  reason  to  expect  that  the  trajectory  of  E^(t) 
is  at  all  times  in  a  classically  allowed  region.  Moreover 
tunneling  is  so  much  associated,  in  our  mind,  with  barrier 
penetration  that  the  above  result  is  somewhat  surprising. 

It  is  interesting  to  note  that  in  Fig.  6f,  corresponding  to 
the  lowest  incident  kinetic  energy,  the  packet  G^  cannot  penetrate 
behind  the  adiabatic  barrier.  The  Ej(t)  curve  comes  down  towards 
the  surface  and  it  is  turned  around  by  it.  It  appears  that 

the  particle  emerges  from  the  surface  in  an  ionic  state,  but  this 
is  not  the  case.  As  shown  in  Fig.  7  the  probability  that  G^  e.xits 
-  for  the  kinetic  energy  used  in  Fig.  6f  -  is  initially  very  small 
and  it  goes  to  zero;  in  contrast  in  the  cases  6a-e,  in  which 
tunneling  takes  place,  the  probability  ( t )  grows  to  finite 
values  (see  Fig.  7). 


III. 4  Two  states  versus  one  state  representation 

Several  problems  in  surface  science  are  conveniently 


represented  in  terms  of  a  two  state  model.  An  obvious  example  is 

alkali  adsorption  on  metals  where  the  adsorbed  alkali  is  ionic  and 

the  incident  one  is  neutral .  The  two  state  character  of  the 

system  is  revealed  by  desoprtion  experiments  in  which  both  ions 

and  neutrals  are  observed.  Another  example  is  dissociative 

adsorption,  which  can  be  thought  of  as  a  transition  from  a  H  -Me 
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state  to  a  2H-Me  state  (Me  means  metal  surface) .  In  this  case 
no  2H  desorption  has  been  observed  to  suggest  that  a  2H-Me  state 
is  involved  in  dynamics.  Therefore  while  one  might  prefer  to  use 
a  two  diabatic  state  description  of  scattering  or  desorption  it 
is  not  at  all  clear  that  a  two  adiabatic  states  description  is 
necessary . 

It  is  therefore  of  interest  to  establish  some  criteria 
whether  one  or  two  state  descriptions  are  necessary.  A  reasonable 
condition  can  be  obtained  by  comparing  the  gap  between  the  two 
lowest  adiabatic  states  with  the  kinetic  energy  of  the  incoming 
particle:  if  the  latter  is  much  smaller  only  one  state  is  needed. 
Of  course,  when  the  two  surfaces  are  unknown,  which  is  generally 
the  case,  one  would  like  to  have  some  observable  effects  which 
signal  the  presence  of  the  second  state.  In  order  to  see  if  such 
effects  exist  we  carried  out  calculations  with  a  Hamiltonian  in 
which  the  adiabatic  surfaces  are  similar  to  those  used  in  the 
previous  section  (Fig.  4),  but  the  gap  between  them  is  smaller 
(Fig.  8).  This  is  obtained  by  using  the  Hamiltonian  II  defined  by 
Eqs .  (  II. 1-3)  and  the  parameters  of  Table  I  with  p=0.1  instead  of 
8=0 . 3 . 

In  Fig.  9  we  show  the  energies  E^(R(t))  and  E2{R(t))  of  the 
ionic  and  neutral  packets,  respectively.  The  incident  kinetic 
energies  are  labelled  by  the  Massey  parameter  X;  the  corresponding 
kinetic  energy  and  the  difference  between  it  and  the  height  of  the 
barrier  are  given  in  Table  II.  The  values  of  X  are  chosen  to  give 
the  same  values  for  the  barrier  height  minus  the  kinetic  energy  as 
those  used  in  the  previous  section. 

The  most  interesting  result  is  that  lowering  the  upper 
adiabatic  state,  to  get  it  closer  to  the  bottom  one,  leads  to 
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higher  neutral  population.  This  happens  as  if  the  upper  state 
helps  reflect  the  Incident  packet.  This  effect  can  be  seen  by 
comparing  the  probabilities  shown  in  Fig.  10  to  those  shown  in 
Fig.  7. 


The  calculations  presented  so  far  have  split  the  incident 
packet  once,  when  it  approached  the  curve  crossing  region.  As  a 
result  the  newly  formed  packet  oscillates  in  the  ionic  state 
forever.  In  reality  the  trapped  packet  splits  off  a  neutral 
packet  whenever  it  approaches  the  curve  crossing  point  and  as  a 
result  the  amplitude  of  the  ionic  packet  dwindles,  and  a 
succession  of  outgoing  neutral  packets  are  created.  This  behavior 
is  the  GWP  description  of  the  "predissociation"  of  the  ionic 
state . 

In  order  to  test  whether  SGWP  generates  such  behavior  we 
carried  out  the  SGWP  calculation  shown  schematically  in  Fig.  11: 
in  Fig.  11(a)  packet  2  is  split  to  generate  packet  1;  after  that 
packet  2  is  turned  around  and  leaves  the  surface  while  1  moves  in 
the  ionic  well  (Fig.  11b);  when  packet  1  approaches  the  crossing 
region  again,  we  use  SGWP  to  generate  a  new  neutral  packet  2'; 
this  leaves  the  surface,  while  1  turns  back  towards  it;  this 
succession  is  repeated  until  the  amplitude  of  1  becomes 
negligible . 

An  idealized  time  of  flight  (TOF)  measurement  applied  to 
this  process  will  give  peaks  corresponding  to  the  arrival  times  of 
the  successive  neutral  packets,  separated  by  a  time  comparable 
with  the  period  of  the  motion  of  G^  in  the  ionic  well.  Such  a 
result  corresponds  to  an  experiment  in  which  the  incident  state  is 
a  packet  that  is  well  localized  in  space.  If  the  experiment  is 
carried  out  with  an  incident  state  which  is  close  to  a  plane  wave 
we  must  describe  it  by  using  a  train  of  several  incident  packets; 
the  succession  of  emerging  neutral  packets  will  then  be  continuous 
in  time.  Moreoever  since  the  time  resolution  of  TOF  detectors  is 
poor,  as  compared  to  the  period  of  the  motion  in  the  ionic  well, 
the  time  "granularity"  given  by  GWP  is  not  a  practical 
shortcoming . 

The  results  of  this  multiple  splitting  calculation  are  shown 
in  Fig.  12.  The  Hamiltonian  is  that  described  in  Fig.  1 ,  X  =  42 
and  p  =  0.3.  The  trajectories  of  the  centers  of  various  packets 


are  plotted  in  Fig.  12a.  We  see  that  the  ionic  packet  is 
formed  at  t=0  and  reaches  the  crossing  point  at  t  ~  6;  later 
turns  around  and  leaves  the  surface,  while  G^  oscillates  in  the 
ionic  well.  As  G^^  moves  away  from  the  surface  a  new  neutral 
packet  G-  is  created  at  t  2s  9,  before  G.  reaches  the  crossing 

2  I  A 

point.  Immediately  G^^  leaves  the  surface  while  G^^  continues  to 
oscillate,  creating  when  it  approaches  the  crossing  point 
again.  The  probability  that  the  particle  is  ionic  is  shown  in 

t  II 

Fig.  12b.  Clearly  the  creation  of  6^  and  G^  diminishes  Pj^(t),  as 
expected.  It  is  not  surprising  that  MTA  give  very  poor  results 
for  this  case. 

We  found  this  procedure  to  be  rather  stable  with  respect  to 
the  precise  point  where  the  packets  are  split.  If  the  splitting 
is  done  too  early  the  amplitude  of  the  new  packet  goes  to  zero;  i 
it  is  done  too  late  the  result  depends  on  the  splitting  point. 
There  is  a  strip  of  about  lA  width,  ahead  of. the  crossing  region, 
in  which  splitting  can  be  done  without  affecting  the  results. 


V. 
V.  1 


THE  EFFECT  OF  DISSIPATION  BY  PHONONS 
Introductory  Remarks 

The  problem  of  curve  crossing  or  tunneling  In  quantum 

systems  coupled  to  a  heat  bath  is  a  topic  of  much  current 
1 6 

interest;  one  would  like  to  formulate  a  theory  which  predicts 

the  manner  in  which  energy  dissipation  and  dephasing  by  thermal 

fluctuations  change  tunneling  rates.  The  present  model  is  also 

relevant  -in  surface  science-  to  the  problem  of  sticking  induced 

by  curve  crossing  followed  by  inelastic  interactions  between  the 

incident  particle  and  the  lattice.  Both  phonons  and  electron  hole 

pair  excitations  may  be  important  and  both  can  be  treated  by  a 

model  in  which  the  incident  particle  is  coupled  strongly  to 

31a 

Independent  bosons . 

In  this  section  we  describe  a  theory  in  which  the  lattice 

motion  is  coupled  to  the  quantum  degrees  of  freedom  through  a  mean 

trajectory  approximation.  This  permits  us  to  combine  the  GWP 

1 4a  5  2 

quantum  dynamics  of  the  particle  with  a  classical  Langevin 

equation  describing  the  motion  of  the  lattice.  Other  models  which 

couple  curve  crossing  dynamics  to  stochastic  variables  simulating 

53 

a  heat  bath  have  been  presented  in  the  literature. 

V . 2  The  model 

We  consider  a  particle  colliding  with  a  one  dimensional  atom 
chain  (Fig.  13)..  The  Hamiltonian  is 


H  =  -(h^/2M)  3^/3R^  -  (Fi^/2m)3^/3r^  -  (h^/2m)  Z3  3q^+H{  r  ,  R ,  x  ) 

(V.  1  ) 

+  V^(r,{q^)) 


The  lattice  atoms  are  divided  here  into  two  groups:  the 
primary  atoms  whose  coordinates  are  denoted  r  and  the  secondary 
atoms  having  the  coordinates  q.  The  concrete  example  used  here 
has  only  one  primary  atom.  The  independence  of  the  electronic 
Hamiltonian  H(r,R,x)  on  the  positions  of  the  secondary  atoms 
implies  that  only  the  atoms  in  the  primary  zone  interact  with  the 
incident  particle.  We  assume  that  the  lattice  is  harmonic  (i.e., 
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the  lattice  potential  energy  V^(r,{q^})  is  quadratic)  but  do  not 
linearize  the  interaction  between  the  incident  particle  and  the 
primary  zone  atoms. 

V . 3  The  Hartree  approximation 

We  make  now  the  Hartree  approximation,  which  assumes  a  wave 
function  of  the  form 

4i  =  Xj^(r,  {q^}  ;t)  CG^(R.t)<t^{R.r.x)  +  ( R ,  t )  ^ ,  x )  ]  (V.2) 

where  are  the  electronic  wave  functions  (1  is  neutral  and  2  is 

ionic)  ,  Gj^  and  are  the  nuclear  wave  functions  of  the  incident 

particle,  and  is  the  wave  function  of  the  lattice. 

Using  the  minimum  error  method^®  and  the  properties  (II. 4) 
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of  the  electronic  wave  functions  leads  to 

ih3G^/3t  =  -(h^/2m)3^G^/3R^  +  {V.3) 


for  i  =1  or  2  and  j  i ,  and  to 

ih3Xj^/3t=-(h^/2m)3^X^/3r^-(h^/2m)Z3^X^/3q^+[Vj^(r,  (q.  )+V^^^(r)  ]X^ 


Here 


(V.4) 


(q^)  )H.  j{R,r)Xj.  (r.  {q.  }  )dr  dq. 


(V.5) 


2  2 

V_..(r)  =  Z  Z  J-G  (R)  H  (R,r)G  (R)dR 
i=l  j=l  i  J 


(V.6) 


H^j(R,r)  =  J'*^(R,r,x)  H(  R ,  r  ,  x )  <t  ^  ( R ,  r ,  x )  dx . 


(V.7) 


Because  we  made  the  Hartree  approximation  the  primary  atom 
interacts  with  the  incoming  particle  through  the  effective 
potential  (V.6).  This  is  the  sum  of  the  ionic  potential  Hj^^(R,r) 


t 


averaged  over  the  distribution  G ^  ( R,  r ;  t )  G (R,  r ;  t )  (which  is 


proportional  to  the  probability  that  the  particle  is  an  ion 
located  at  R)  plus  the  neutral  potential  averaged  over 


the  distribution  G ( R,  r  ;  t )  G ^^  ( R,  r ;  t )  (which  is  proportional  with 


the  probability  that  the  particle  is  neutral  and  is  located  at  R] 

« 


plus  the  interference  term  2ReXdRGj^(R,  r ;  t )  G2(R,r;t) 


To  do  better  than  the  Hartree  approximation  we  would  have  to 
use  the  wave  function 


i)((R,r{q}  ,x;t) 


=  X^({q.)) 


2 

1 


1=1 


G^  ^(r,R;t)*^(R,r,x) 


(V.8) 


which  Incorporates  the  correlated  quantal  motion  of  the  primary 
zone  atoms  and  that  of  the  incident  particle;  the  Hartree 
approximation  is  made  for  the  secondary  lattice  atoms  only. 

It  is  Instructive  to  think  under  what  conditions  we  can  hope 
that  the  Hartree  approximation  is  satisfactory.  This 
approximation  simplifies  the  force  exerted  by  the  incident  atom  on 
the  primary  atom,  by  averaging  it  over  the  state  of  the  incident 
particle.  If  the  particle  state  Z  G^  is  predominantly  ionic 
than  the  force  exerted  on  the  primary  atom  resembles  that  caused 
by  the  ion;  in  the  opposite  case  the  interaction  resembles  that 
exerted  by  the  neutral .  In  the  intermediate  case  the  force  is  the 
average  of  the  ionic  and  neutral  forces,  with  the  corresponding 
quantum  weights.  For  the  present  problem  one  may  assume  that  the 
lattice-atom  energy  transfer  and  dephasing  are  most  efficient  at 
the  moment  of  the  impact  and  that  the  repulsive  part  of  the 
potential  is  most  important.  Since  the  repulsive  walls  for  ion 
and  neutral  are  similar  it  may  not  be  essential  that  the  effective 
potential  mixes  them  in  exactly  the  right  way.  On  the  other  hand, 
if  electron  hole  pair  excitations  and  the  polarization  of  the 
electron-gas  (i.e.  "image  effect")  are  important  we  might  want  to 
use  the  non-Hartree  wave  function  Eq  (V.8).  This  will  treat  the 
coupling  of  the  ion  to  the  electron  gas  on  a  different  footing 
than  the  coupling  of  the  neutral. 
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V. 4  The  mean  trajectory  approximation  for  lattice  motion 

The  equations  (V.3-6)  can  be  further  simplified  by  assuming 
that:  (a)  the  lattice  wave  function  is,  throughout  the  collision, 

Gaussian  in  the  variables  q  and  r;  and  that  (b)  the  mean  square 
displacement  of  the  variable  r  in  this  Gaussian  state  is  very 
small.  Assumption  (b)  allows  us  to  perform  the  integral  in  Eq. 
(V.5)  by  the  steepest  descent  method,  to  obtain 

=  H^j(R,r{t))  ,  (V.9) 

where  r(t)  is  the  mean  position  of  the  primary  atom  at  time  t: 

M 

r(t)  =  Xdr  II  dq  X-(r,{q):t)  rX^ ( r , {q} ; t ) .  (V.IO) 

i=l  1  ^  ^ 

Due  to  Eq.  (V.9)  and  (V.3)  the  wave  functions  G^  and  G^  no 
longer  depend  on  the  lattice  wave  function,  but  only  on  the  mean 

e 

primary  atom  position  r(t).  Prom  Heller's  work  we  know  that  if 
the  Hamiltonian  is  harmonic  in  the  variables  {q.},  and  the 
Gaussian  wave  function  X^  is  narrow  with  respect  to  r,  then  it 
follows  that  r(t)  and  q^(t)  satisfy  classical  equations  of  motion: 

mr(t)  =  -3V^,^^/3r  { t)  -  mo^  (  r  ( t ) -q^  ( t )  )  ,  ^  (V.ll) 

w“^mq^(t)  =  q^_^(t)-2q.(t)  +qi  +  i(t),  1  =  1,2 -  (V.12) 

with 

q^lt)  =  r{t) 

The  mean  trajectory  approximation  decouples  thus  the  wave 
functions  G^,i  =  l,2  from  the  wave  function  X^^ ;  the  only  lattice 
information  that  is  needed  is  the  time  evolution  of  the  center  of 
the  packet  representing  the  primary  lattice  atom;  the  latter  moves 
classically  on  a  mean  potential  which  is  the  average  interaction 
between  the  incident  particle  and  the  primary  atoms  over  the 


particle  wave  function. 

To  assess  the  validity  of  the  mean  trajectory  approximation 
we  must  quantify  the  demand  that  the  Gaussian  is  narrow.  Since 
this  assumption  is  used  to  justify  the  integration  of  (V-5)  by  the 

Laplace  method,  the  width  l(t)  of  the  Gaussian 

« 

g(r)=JXj^(r,qj.  .qj^)X^(r,q^.  .qj^)  must  be  smaller  than  the 

length  L  over  which  the  function  H^j(r,R)  changes  with  r. 
Furthermore,  for  r(t)  to  move  classically  under  the  influence  of 
the  potential  must  have^*^'^^ 

l^(t)<<4[3Vgj^(r(t) )/3r(t) ]"^[3^V^^^(r(t) )/3r(t)^]  =  n ( t ) .  An 

order  of  magnitude  estimate  of  the  width  l(t)  of  the  Gaussian  g(r) 

2 

is  the  mean  displacement  of  the  primary  atom  1  ''(h/mu) 
coth[h<i)/2kgT]  . 

As  a  side  remark  we  note  that  this  expression  for  1^  does 

not  lead  to  a  divergence,  in  the  present  theory,  for  w-»0 .  Indeed 
2 

as  u-*0  and  1  -*«  the  Gaussian  g(r)  tends  to  I6(r-r(t)).  Since  in 

the  present  theory  the  wave  function  appears  everywhere  in 

integrals  over  frequency  the  quantity  l(w)  is  multiplied  with  the 

density  of  states,  leading  thus  to  expressions  containing 

2  2 

/dwp («)!(«).  Since  p(«)  =  3u  (for  a  Debye  model)  l(w)p(w)  is 

proportional  to  u  as  «-»0.  Thus  in  the  w-^O  limit  the  integrals 
through  which  l(w)  appears  in  the  theory  are  well  behaved. 

We  can  summarize  the  conditions^under  which  the  Gaussian 
g(r)  is  sufficiently  narrow  by  requiring 

(h/mu)coth(hu/2kgT)<<max[L(t)^,il(t)^]  =  . 

In  the  high  temperature  limit  hu/k„T<<l  and  coth ( hw/ 2kT ) -►2k„T  hw . 

2  2°  ° 

Thus  the  condition  ( 2h  /mX  )<<k_T  ensures  the  validity  of  the  mean 

o 

trajectory  approximation  in  the  high  temperature  limit.  This 

result  makes  sense,  since  in  the  high  temperature  limit  the 

important  frequency  in  the  system  is  k_T/h.  We  can  use  it  to  form 

2  “ 

the  action  A  a  mX  k^T/h  by  using  X  as  a  typical  length  of  the 
potential.  Then  a  classical  trajectory  limit  should  be  valid  if 
h/A<<l,  which  leads  to  the  high  temperature  condition  obtained 
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above . 

V . 5  A  summary  of  the  equations 

After  we  made  the  mean  trajectory  approximation  for  the 
motion  of  the  primary  atom  and  took  the  classical  limit  for  the 


motion  of  all  the  lattice  atoms,  the  nuclear  wave  functions 
satisfy 


lh3G^(R, t)/3t=-(h^/2M)3^G^(R,t)/3R^+H^^(R.r( t) )G^(R,t) 

+  H^j(R,r(t)  )Gj(R,t)  ,  i,j=1.2  and  i?ij.{V.13) 

These  equations  are  solved  by  the  MEM-SGWP  method. 

To  compute  r(t)  we  use  the  Eqs.  (V. 11-13)  and  the  Lange vin 
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equation  method  proposed  by  Adelman  and  Doll.  The  version 
implemented  here  is  the  ghost  atom  method  developed  by  Garrison 

e  e  e  C 

and  Adelman  and  Shugard,  Tully.and  Nitzan.  This  solves 
mr(t)  =  -3V^^^/3r(t)-moj^(r-r^)  (V.14) 

and 

••  2  2 

m^glt)  =  -m«  (rg-r)-m«grg  -  mydr^/dt+F ( t )  (V.15) 


In  choosing  the  friction  coefficient  y  and  the  ghost  atom 
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frequency  we  follow  Garrison  and  Adelman.  The  fluctuating 
force  has  the  correlation  function 


<F(t)F(0)>  =  2mkgTy6(t) 


(V. 16) 


which  assumes  the  high  temperature  (i.e.  classical)  limit  when  the 
thermal  average  is  performed.  Some  quantum  effects  can  be  put 
back  into  the  theory  by  replacing  2k_T  with  hwcoth( hw/2knT) .  The 

D  D 

qualitative  effect  of  using  the  expression  (V.16)  is  that  at  low 
temperature  (i.e.  hu/k^T  >1)  it  produces  more  noise  than  the  rea.T. 
quantum  system  does.  Furthermore,  Eq.  (V.16)  does  not  give 
different  weights  to  energy  loss  (i.e.  Stokes)  and  energy  gain 
(i.e.  anti-Stokes)  processes.  This  should  introduce  very  large 


errors  at  the  lowest  temperatures  where  the  Stokes  processes  have 

high  intensity  and  that  of  anti-Stokes  processes  is  almost  zero. 

As  Tully's  extensive  work  has  demonstrated  the  method  based  on  Eq . 

(V.16)  should  be  reliable  at  high  temperatures  (i.e.  hu<<k_T) . 

o 

V.  6  Illustrative  numerical  results 

We  present  results  of  several  calculations  carried  out  by 
solving  Eqs .  (V. 13-16).  The  lattice  parameters  are 

_3 

Wjj=l. 43x10  a.u.  and  the  mass  of  Ni .  The  projectile  mass  is  that 

of  Na  and  the  matrix  (R,r(t))  is  defined  by  Eqs.  (II. 1-3)  with 

R  replaced  with  R-r(t);  the  values  of  the  parameters  are  those 

listed  in  Table  I  under  Hamiltonian  I.  The  coupling  strength  is 

^=0.3  and  the  incident  kinetic  energy  is  given  by  X=43  (the 

2  -2 

kinetic  energy  is  =  89. 9X  eV) . 

In  Figs.  13a  we  show  the  trajectory  followed  by  the  packets 
(one  crossing  only)  in  the  case  when  the  lattice  is  rigid.  The 
effect  of  allowing  the  lattice  to  move,  at  T=0K  is  shown  in  Fig. 
13b,  where  the  amplitude  of  the  ionic  packet  decays  as  a  result  of 
phonon  excitations.  The  zero  point  energy  of  the  lattice  was 
included  in  the  classical  Langevin  equation.  In  Fig.  13c  we  show 
the  energy  Ej^(t)  of  the  ionic  packet  (see  Eq.  (III. 4))  as  a 
function  of  the  packet  position.  The  arrows  indicate  the  flow  of 
time.  The  behavior  of  the  system  at  non-zero  temperatures  is 
shown  in  Figs.  14a  and  b.  At  lower  temperature  (1=0^^,  where  0j^  is 
Debye  temperature)  the  ion  is  quickly  equilibrated,  and  phonon 
absorption  is  seen;  however  the  amplitude  does  not  decay  as 
continuously  as  in  Fig.  13a  at  T=0 .  For  T=50jj  (Fig.  14b)  the  ion 
has  not  settled  to  a  low  amplitude  even  after  three  oscillations 
in  the  well. 

These  calculations  illustrate  the  fact  that  the  addition  of 
a  thermal  dissipative  channel  drives  the  Gaussian  wave  function  in 
a  manner  that  resembles  classical  mechanics.  The  computation  is 
rather  inexpensive  and  we  estimate  that  it  is  possible  for  three 
dimensional  models. 

VI.  Conclusions 

We  have  presented  a  Guassian  wave  packet  theory,  for  a  curve 
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crossing  system,  which  has  one  packet  per  electronic  state.  The 
most  interesting  feature  of  the  theory  is  that  it  permits 
tunneling,  in  the  sense  that  a  packet  incident  in  a  barrier  splits 
into  one  packet  reflected  by  the  barrier  and  one  which  penetrates 
behind  it.  The  procedure  for  multiplying  the  number  of  packets  in 
the  course  of  propagation  is  general  and  should  be  useful  in 
numerous  other  applications  of  the  GWP  method  to  phenomena  in 
which  wave  function  splits  spontaneously,  in  the  course  of  the 
collision,  into  spatially  disjoint  pieces. 

Another  new  development  is  the  discussion  of  various 
approximations  Involved  in  coupling  the  curve  crossing  system 
whose  properties  are  described  within  a  GWP  formulation  to  a  many- 
body  classical  system,  such  as  a  lattice.  The  use  of  Hartree 
approximation  permits  an  approximate  but  self-consistent  treatment 
of  this  problem.  The  classical  limit  for  the  motion  of  the  many- 
body  system  is  not  necessary;  the  use  of  the  GWP  method  to 
propagate  the  bath  variables  quantum  mechanically  is  feasible,  as 

e  Q 

demonstrated  recently  by  Singer  and  Smith. 

The  numerical  studies  presented  here  explored  the  stability 
of  the  method  and  studied  the  qualitative  behavior  of  the  packets. 
We  are  now  developing  an  exact  method  for  solving  this  problem  and 
hope  to  report  soon  a  comparison  between  the  results  obtained  with 
SGWP  and  "the  exact  ones.  If  the  SGWP  method  turns  out  to  have 
satisfactory  accuracy  than  it  can  be  used  for  three-dimensional 
studies  for  which  exact  calculations  are  much  more  difficult. 

We  had  two  reasons  for  undertaking  this  study.  First,  we 
were  curious  about  the  behavior  of  the  packets  in  a  strongly  non- 
classical  situation,  such  as  a  two  state  system  at  energies  for 
which  tunneling  plays  an  important  role.  Second,  we  are  searching 
for  a  method  of  doing  very  inexpensive  three-dimensional  curve 
crossing  calculations  for  a  quantum  system  imbedded  in  a  classical 
heat  bath;  since  the  behavior  of  the  bath  must  be  generated  by 
repeated  use  of  classical  molecular  dynamics  the  ability  to 
calculate  the  quantum  part  cheaply  {but  reliably)  is  of  paramount 
importance.  The  GWP  method  seems  a  natural  candidate  for  this 
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APPENDIX 

The  Equations  used  for  the  MEM-SGWP  propagation  of  the 
57 

packets  are: 


h  “  ^1^1  =  ~  Pl/2m  +{[V^.(2)+V^j(2)]  M^(2)  - 

[V^^(0)+V^j(0) ]M^(4) 

+  2a^/m  =  B . {M^ ( 2 ) [ V . ^ ( 0 )  +  V^j(O)] 

-  M^(0)  [V^^(2)  +  V^j(2)]} 

^1  "  -V^^(1)M^(2)“^  +  (M. (2)lma^)"^  Im[a*V^  (1)] 


-1 


=  P^/m  +  {2Ima^  M^(2)}  “  Im  V^^d) 


where 

M^(n)  =  ;G^(R)*(R-R^(t) )"g^(R)  , 


V.jin)  =  ;G^(R)*(R-R^(t)  )^^j(R)Gj(R)dR 

B^  =  {M(4)M^(0)  -  M(2)^}"^ 
and  1=1,2,  j=i,2,  lx j . 


Also  note  that 


(A.l) 


(A. 2) 

(A. 3) 

(A. 4) 

(A. 5) 

(A. 6) 
(A. 7) 


. 


Nv 

■,v\' 

V-N 


i;A 

v"/ 

'A 

•  *.'  1 

d 


-1 


-  V^^(l)M^{2)  =  -XdRG^G^( 3H.^/3R)/XdRG.G. 


gf-  {XdRG*H..  G./XdRG*G.} 


(A. 8) 


r  j 


and  that 


B  (M  (2)V  (0)-V  (2)M  (0)  }=  -  (1/2)  {XdG*H  G  /;dRG*G  } 

3R^ 

(A. 9) 


f 


■*  ■  ^  •-  'V  •-  •-  '  ■  ^  .^,1  ^,1  — 


Table  I.  The  values  of  the  parameters  used  In  Eqs.  (11.1-3) 

to  define  the  Hamiltonians. 


Hamiltonian  I  Hamiltonian  II 


atomic  units 

other 

atomic  units 

other 

u 

0.025H^®^ 

0.68  eV 

0.184H^®^ 

5.0  eV 

o 

Ae 

O.OOSH^^^ 

0.136  eV 

0.147H^®^ 

4.0  eV 

2.5  a<^) 

1.32  A 

2.5  a^^^ 

1 .32A 

Ro 

5.0  a<^) 

2.64  A 

5.0  a<^) 

2.64A 

R 

12.5  a^^^ 

6.6  A 

9.0  a<^> 

4.75A 

c 

p 

0.1  or  0.3 

0.1  or  0.3 

0.1  or  0.3 

0.1  or  0.3 

m 

42300^®^ 

22.991^*^^ 

12 , 800a. u. 

(d)  e  g4(d) 

(a) 

(b) 

(c) 

(d) 

Hartree 

a  is  Bohr  radius 

the  mass  of  Na 

the  mass  of  Li 

r 


Table  II. 


The  energies  of  the  Incident  packets.  For  the 
definition  of  column  headings  see  the  text. 


Incident  energy-barrier  height 


Pg/Zm 

LP 

g/Sm 

8=0.3 

8=0.1 

90 

2.91 

1.39 

X 

10-2 

— 

9.21  X 

10-^ 

100 

2.35 

1.76 

X 

10-2 

7.60  X  10-^ 

3.70  X 

10-1 

110 

1.95 

1.42 

X 

10-2 

3.48  X  10-^ 

-4.16  X 

10-1 

120 

1.63 

1 . 17 

X 

10-2 

3.50  X  10-2 

-3.54  X 

10-1 

130 

1.39 

9 . 82 

X 

10-3 

-2.10  X  10-^ 

-5.99  X 

10-1 

140 

1.20 

8.34 

X 

10-3 

-4.05  X  10“^ 

-9.49  X 

10-1 

150 

1.05 

7.25 

X 

10-3 

-5.60  X  10-^ 

— 

— 

eV 

eV 

eV 

eV 

(a)  Massey  parameter  (Section  II, 4) 


Table  III. 


The  probability  that  a  particle  starting  in  the 
state  2  emerges  in  state  1  (one  crossing  only) . 


p  =  0.3 


X 

LHA- 

SGWP 

MEM- 

SGWP 

MTA 

40 

0.954 

0.655 

1.0 

35 

41 

0.979 

0.588 

1.0 

36 

42 

0.954 

0.500 

1.0 

37 

43 

0.887 

0 . 404 

1.0* 

38 

44 

0.750 

0.308 

1.0* 

39 

45 

0.622 

0.219 

1.0* 

40 

+ 

The  packets  turn 

around 

before 

region. 


P  = 

0.1 

LHA- 

MEM- 

MTA 

SGWP 

SGWP 

0.429 

0.404 

0.807 

0.409 

0.382 

0 . 888 

0 . 372 

0.334 

0.978 

0.341 

0.286 

+ 

1 

0.370 

0.235 

+ 

1 

0.395 

0 . 180 

t 

1 

reaching  the  curve  crossing 


FIGURE  CAPTIONS 


Fig.  1. 


Fig.  2. 


Fig.  3. 


Fig.  4. 


The  diabetic  potentials  i  =  1,2  (see  Eqs. 

II.  2),  the  coupling  H^j(R)  (dash-dotted  line,  right 
hand  side  scale)  and  the  adiabatic  potentials  H®^(R), 
(dashed  lines)  i=l,2.  The  parameters  are  those  for 
Hamiltonian  I  in  Table  I,  with  p=0.3.  The  adiabatic 
and  diabatic  potentials  overlap  everywhere  except  in 
the  crossing  region. 

The  trajectories  of  the  centers  R^  (dashed  line)  and 

R^  (full  line)  of  the  ionic  and  neutral  packets, 

respectively,  obtained  by  MEM-SGWP.  The  dash-dotted 

line  is  the  MTA  result.  The  Hamiltonian  I  with  p=0,3 

(see  Table  I  for  parameter  values  and  Fig.  1  for 

plots)  was  used  for  all  calculations.  The  Massey 

parameters  for  the  incident  G2  packets  are:  (a) 

X=37 ;  (b)  X=38;  (c)  X=42  and  (d)  X=43.  The  kinetic 
2  -2 

energy  is  mVQ/2  =  3.30  X  Hartree  = 

89.9  X"^  eV. 

The  time  evolution  of  the  probability  that  a  neutral 
particle  approaching  the  wall  remains  neutral.  The 
full  lines  are  the  results  obtained  with  MEM-SGWP  and 
the  dashed  lines  those  given  by  MTA.  The  parameters 
are  those  used  in  Fig.  2. 

The  adiabatic  potentials  H^j^(R)  and  (full 

lines)  and  the  diabatic  coupling  (dashed  - 

dotted  line),  given  by  Eqs.  (II. 1-3)  and  the 
parameters  for  Hamiltonian  II  of  Table  I,  with  3=0.3. 
The  diabatic  curves  H^j^  and  H^j  are  the  dashed  lines. 
They  differ  from  the  diabatic  curves  only  in  the 
crossing  region.  The  arrows  on  the  right  hand  side 
indicate  the  kinetic  energy  of  the  incident  packets. 


.■» '■T'  iT" I 


Fig.  5, 


Fig.  6, 


Fig.  7 


Fig.  8. 


For  conversion  from  the  Massey  parameter  X  to  the 
kinetic  energies  see  Table  II, 

The  trajectories  Rj(t)  and  RgCt)  of  the  centers  of 
the  ionic  (full  lines)  and  neutral  (dashed  lines) 
packets  and  respectively.  The  calculations 

are  for  the  Hamiltonian  II  defined  by  Eqs .  (II. 1-3) 
and  the  parameters  given  in  Table  I,  with  p=0.3.  The 
incident  Massey  parameters  are:  (a)  X=100;  (b)  X=110; 

(C)  X=120;  (d)  X=130;  (e)  X=140;  (f)  X=150.  The 

calculations  (a)  -  (c)  correspond  to  kinetic  energies 
above  the  barrier,  the  others  are  below.  For 
conversion  of  X  to  kinetic  energies  see  Table  II. 

The  adiabatic  energy  surfaces  and  (dash- 

dotted  lines)  and  the  packet  energies  Ej(t)  and  E2 ( t ) 
defined  by  Eq.  (III.4)).  The  calculations  are  for 
the  Hamiltonian  II  defined  by  Eqs.  (II.  1-3)  and  the 
parameters  given  in  Table  I  with  p=0.3.  The  distance 
R  is  the  position  of  the  center  of  the  corresponding 
packet.  The  arrows  on  the  curves  indicate  the  flow 
of  time.  Various  curves  correspond  to  (a)  X=100;  (b) 

X=110;  (c)  x=120:  (d)  X=130;  (e)  X=140;  X=150. 

The  probability  that  a  particle  starting  in  the 
neutral  state  2  ends  up  in  the  state  1 .  Various 
curves  correspond  to  different  Massey  parameters  \. 
The  correspondence  between  X  and  the  incident  kinetic 
energy  is  given  in  Table  II. 


The  matrix  elements  H,.(R)  (full  lines)  and 

ii 

H^ j ( dashed-dotted  lines)  defined  by  Eqs.  (II. 1-3)  and 
the  parameters  for  Hamiltonian  II  given  in  Table  I. 
The  value  of  p  is  0.1.  The  adiabatic  energy  curves 
H®  are  also  given.  The  incident  kinetic  energies 


used  in  subsequent  calculations,  labeled  by  the 
values  of  X  (see  Table  II  for  conversion)  are  marked 
on  the  graph. 

Pig.  9.  The  adiabatic  energy  surfaces  of  Fig.  8  (dashed- 

dotted  lines)  and  the  energies  E^(R(t))  and 
of  the  ionic  and  neutral  packets,  respectively  (see 
Eq.  (III. 4)  for  definition),  (a)  X=90;  (b)  X=100;  (c) 

X=110;  (d)  X=120;  (e)  X=130;  (f)  X=140.  The  Massey 

parameter  X  indicates  the  incident  kinetic  energy 
( see  Table  II )  . 

Fig.  10.  The  probability  Pj(t)  that  a  neutral  particle  is 

ionized,  for  the  Hamiltonian  shown  in  Fig.  8  the 
curves  corresponding  to  different  incident  kinetic 
energies  are  labelled  by  X.  For  the  relationship 
between  X  and  kinetic  energy  see  Table  II. 

Fig.  11.  Schematlcal  representation  of  successive  splitting  of 

packets.  See  the  text  for  explanation. 

Fig.  12.  The  behavior  of  the  Gaussian  packets  when  multiple 

splitting  is  performed.  The  calculation  is  done  with 
Hamiltonian  I  (see  Eqs .  (II. 1-3)  and  Table  I),  p=0.3 
and  X=42  (mv^/2  =  89.9  X  ^).  (a)  The  evolution  of 

the  centers  of  the  neutral  packets  ,  R^  -  Rj 
that  of  the  ionic  packet  R^ .  (b)  the  probability 

Pj(t)  that  the  system  is  ionic. 

Fig.  13.  Schematic  representation  of  the  lattice  atom  system 

and  its  parameters. 

Fig.  14.  The  results  of  a  calculation  using  the  Hamiltonian  I 

(Table  I)  X=43  (mv^/2  -  89.9  x"^  eV)  and  p=0 . 3 , 
coupled  through  a  mean  trajectory  Langevin  equation 


to  a  Ni  lattice  with  (i>j^=l .  43x10  a.u.  (a)  The  motion 
of  the  centers  of  the  ionic  (R^)  and  neutral  (R^) 
packets  for  a  rigid  lattice;  (b)  the  same  for  a 
lattice  at  T=0K  (c)  the  energy  of  the  above  packets 
at  T=0K. 

Same  as  Fig.  13b  except  that;  (a)  7=0^^ :  and  (b) 

T=50_,  where  0_  is  the  Debye  temperature. 
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